function res = covadj (x, Ga, Gc)

   ## usage: res = covadj (x, Ga, Gc)
   ##
   ## Adjust linear system to remove model covariance bias between control (Gc) and ana (Ga).

   global PAR

   f = dbstack ;

   J = sum(!isnan(x)) >= PAR.Ccases ; x = x(:,J) ;

   if nargin < 2

      c = cov(x) ;

      [g p] = chol(c) ;
      if (p != 0)
	 [g p] = chol(bend(c)) ;
	 if (p != 0)
	    save /tmp/tt.ob c
	    error("covadj: matrix not pd, saving to /tmp/tt.ob\n") ;
	 endif 
      endif

      res = g ;      

   else

      C = Gc \ Ga ;
      printf("%s: size of covariance adjustment: %3.1f\n", f(1).name, norm(C)) ;
      res = nanmult(x, C) ;

   endif

endfunction
